/*
 * =====================================================================================
 *
 *       Filename:  logistic.h
 *
 *    Description:  logistic regression
 *
 *        Version:  1.0
 *        Created:  2009年06月26日 09时31分48秒
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Ying Wang (WY), ywang@nlpr.ia.ac.cn
 *        Company:  Institute of Automation, Chinese Academy of Sciences
 *
 * =====================================================================================
 */
#include "ncmatrix.h"
#include "ncvector.h"
#include "mathutils.h"
#include "blas.h"
/**
 * w=[w0,w1,...,wd], p(C1|x) = \frac{1}{1+\exp(-w^Tx)}, x=[1,x1,...,xd]
 *
 */
NCvector<double> logistic_regression(const NCmatrix<double> &data1, const NCmatrix<double> &data2)
{
	int n,d,dim = data1.column(), N1= data1.row(), N2= data2.row();
	double tmp;
	NCmatrix<double> data(N1+N2,dim+1);
	NCvector<double> t(N1+N2);
	NCvector<double> w(dim+1,1.);
	NCvector<double> y(N1+N2);
	NCvector<double> grad(dim+1);
	NCmatrix<double> H(dim+1,dim+1);
	NCmatrix<double> RQ(N1+N2,dim+1);
	//NCvector<double> 

	//
	for( n=0; n<N1; n++ )
	{
		data[n][0] = 1.;
		t[n] = 0.;
	}
	for( n=0; n<N2; n++ )
	{
		data[N1+n][0] = 1.;
		t[N1+n] = 1.;
	}
	//
	for( d=1; d<dim+1; d++ )
	{
		for( n=0; n<N1; n++ )
			data[n][d] = data1[n][d-1];
		for( n=0; n<N2; n++ )
			data[N1+n][d] = data2[n][d-1];
	}

	//y
	for( n=0; n<N1+N2; n++ )
	{
		for(tmp=0., d=0; d<dim+1; d++ )
		{
			tmp += data[n][d]*w[d];
		}
		y[n] = sigmoid(tmp);
	}
	//gradient
	grad = data.transpose() * (y-t);
	for( n=0; n<N1+N2; n++ )
	{
		for( d=0; d<dim+1; d++ )
		{
			RQ[n][d] = (1-y[n])*y[n]*data[n][d];
		}
	}

	// H = data.transpose() * RQ w= w- inv(H)*grad
	w  = w - solve(data.transpose() * RQ, grad);

	
	return w;
	//Hessian matrix
}

